import java.util.Random;
import java.util.Arrays;
import java.lang.Math;
import java.io.*;

/*
\item Comparing error (absolute) of each sketch as memory percent changes
\item Varying size of data
\item Plotting actual values as some parameter (location/shape) is varied around actual value
*/

public class OneSampleKSTest
{
    public static void main(String[] args) throws Exception // makes code simpler to avoid try/catch
    {
	if (args.length < 3)
	{
	    System.out.println("Insufficient arguments");
	    System.exit(0);
	}

	OneSampleKSTest ks = new OneSampleKSTest();

	if (args[0].equalsIgnoreCase("M"))  // vary memory
	    ks.varyMemory(args[1], args[2]);
	else if (args[0].equalsIgnoreCase("S"))  // vary data size
            ks.varySize(args[1], args[2]);
	else if (args[0].equalsIgnoreCase("D"))  // vary distance
            ks.varyDistance(args[1], args[2]);
    }

    // cdf function being compared with for varyMemory
    public double cdf(double x)
    {
	// return cdfNormal(x, 0, 1);   // N(0, 1)
	return x;  // U(0, 1)
	//return 1.0 - Math.pow(1.0/x, 1.0);  // pareto dist with scale 1.0, and shape 1.0
    }

    public void varyMemory(String d1, String results) throws Exception // makes code simpler to avoid try/catch
    {
	int n = 10000;
	int datarepeats = 10;
	int repeats = 1;
	Random r = new Random();
	double[] frac = {0.001, 0.0031622, 0.01, 0.031622, 0.1}; //{0.05, 0.1, 0.15, 0.20, 0.25};
	int fracs = frac.length;
	double[][] relErr = new double[4][fracs];
	double[][] absErr = new double[4][fracs];

	for (int i = 0; i < fracs; ++i)
            for (int j = 0; j < 4; ++j)
                relErr[j][i] = absErr[j][i] = 0.0;

	BufferedWriter w = new BufferedWriter(new FileWriter(new File(results)));

	for (int k = 0; k < datarepeats; ++k)
	{
	    // read data from file
            BufferedReader r1 = new BufferedReader(new FileReader(d1 + k + ".dat"));
            n = Integer.parseInt(r1.readLine());
            double[] data = new double[n];
            for (int i = 0; i < n; i++)
                data[i] = Double.parseDouble(r1.readLine());
            r1.close();
	    
	    // created sorted version of the data
	    double[] dataSorted = new double[n];
	    for (int i = 0; i < n; i++)
		dataSorted[i] = data[i];
	    Arrays.sort(dataSorted);
	    
	    // compute the real KS statistic
	    double D_real = 0.0;
	    for (int i = 0; i < n; i++)
	    {
		D_real = Math.max(D_real, (i + 1.0)/n - cdf(dataSorted[i]));
		D_real = Math.max(D_real, cdf(dataSorted[i]) - ((double)i)/n);
	    }
	    System.out.println("\nReal: D = " +  D_real);
	    w.write("\nReal: D = " +  D_real + "\n");
	    
	    int count = 0;
	    for (int l = 0; l < fracs; ++l)
	    {
		double f = frac[l];
		System.out.println("\nf = " + f);
		w.write("\nf = " + f + "\n");
		
		/** test the GK-data structure **/
		GreenwaldKhannaQuantileSummary s = new GreenwaldKhannaQuantileSummary((int)(f * n));
		for (int i = 0; i < n; ++i)
		    s.insert(data[i]);
		// compute the estimated KS statistic from GK sketch
		double D_GK = 0.0;
		for (int i = 0; i < n; i++)
		{
		    D_GK = Math.max(D_GK, (i + 1.0)/n - cdf(s.getQuantile(((double)i)/n)));
		    D_GK = Math.max(D_GK, cdf(s.getQuantile(((double)i)/n)) - ((double)i)/n);
		}
		System.out.print("GK: " + D_GK);
		System.out.print("; " + Math.abs(D_GK - D_real));
		System.out.println("; " + relError(D_GK, D_real));
		w.write("GK: " + D_GK + "; " + Math.abs(D_GK - D_real) + "; " + relError(D_GK, D_real) + "\n");
		relErr[0][l] += relError(D_GK, D_real)/datarepeats;
		absErr[0][l] += Math.abs(D_GK - D_real)/datarepeats;
	    
		/** test the qDigest data structure **/                 
		QDigest s2 = new QDigest((int)(f * n)/3);  // compression factor = space/3                 
		for (int i = 0; i < n; ++i)  
		    s2.insert(data[i]);  
		
		// compute the estimated KS statistic from qDigest sketch      
		double D_QD = 0.0;     
		for (int i = 0; i < n; i++)     
		{ 
		    //System.out.println(((double)i)/n + " " + s2.getQuantile(((double)i)/n));
		    D_QD = Math.max(D_QD, (i + 1.0)/n - cdf(s2.getQuantile(((double)i)/n)));                     
		    D_QD = Math.max(D_QD, cdf(s2.getQuantile(((double)i)/n)) - ((double)i)/n);                 
		}
		System.out.print("QD: " + D_QD); 
		System.out.print("; " + Math.abs(D_QD - D_real));  
		System.out.println("; " + relError(D_QD, D_real));   
		w.write("QD: " + D_QD + "; " + Math.abs(D_QD - D_real) + "; " + relError(D_QD, D_real) + "\n");
		relErr[1][l] += relError(D_QD, D_real)/datarepeats;     
		absErr[1][l] += Math.abs(D_QD - D_real)/datarepeats; 
		
		
		/** test the reservoir sampling data structure **/      
		ReservoirSampler s3 = new ReservoirSampler((int)(f * n));  
		for (int i = 0; i < n; ++i)  
		    s3.insert(data[i]);                
		// compute the estimated KS statistic from reservoir sampling sketch   
		double D_RS = 0.0;    
		for (int i = 0; i < n; i++) 
		{   
		    D_RS = Math.max(D_RS, (i + 1.0)/n - cdf(s3.getQuantile(((double)i)/n)));
		    D_RS = Math.max(D_RS, cdf(s3.getQuantile(((double)i)/n)) - ((double)i)/n);
		}
		System.out.print("RS: " + D_RS);
		System.out.print("; " + Math.abs(D_RS - D_real));  
		System.out.println("; " + relError(D_RS, D_real));  
		w.write("RS: " + D_RS + "; " + Math.abs(D_RS - D_real) + "; " + relError(D_RS, D_real) + "\n");
		relErr[2][l] += relError(D_RS, D_real)/datarepeats;
		absErr[2][l] += Math.abs(D_RS - D_real)/datarepeats; 
	    
		/** Compute value from exact KS on sampled data **/
                double[] sample = s3.getSamples();
		Arrays.sort(sample);
		// compute the estimated KS statistic from sampling
                double D_SA = 0.0;
                for (int i = 0; i < sample.length; i++)
		{
		    D_SA = Math.max(D_SA, (i + 1.0)/sample.length - cdf(sample[i]));
		    D_SA = Math.max(D_SA, cdf(sample[i]) - ((double)i)/sample.length);
		}
                System.out.print("SA: " + D_SA);
                System.out.print("; " + Math.abs(D_SA - D_real));
                System.out.println("; " + relError(D_SA, D_real));
                w.write("SA: " + D_SA + "; " + Math.abs(D_SA - D_real) + "; " + relError(D_SA, D_real) + "\n");
                relErr[3][l] += relError(D_SA, D_real)/datarepeats;
                absErr[3][l] += Math.abs(D_SA - D_real)/datarepeats;
	    }
	}

	for (int i = 0; i < fracs; ++i)
	{
	    System.out.println("f = " + frac[i]);
	    System.out.println("GK: " + relErr[0][i] + ", " + absErr[0][i]);
	    System.out.println("QD: " + relErr[1][i] + ", " + absErr[1][i]);
	    System.out.println("RS: " + relErr[2][i] + ", " + absErr[2][i]);
	    System.out.println("SA: " + relErr[3][i] + ", " + absErr[3][i]);

	    w.write("f = " + frac[i] + "\n");
	    w.write("GK rel: " + relErr[0][i] + "\n");
	    w.write("GK abs: " + absErr[0][i] + "\n");
	    w.write("QD rel: " + relErr[1][i] + "\n");
	    w.write("QD abs: " + absErr[1][i] + "\n");
	    w.write("RS rel: " + relErr[2][i] + "\n");
	    w.write("RS abs: " + absErr[2][i] + "\n");
	    w.write("SA rel: " + relErr[3][i] + "\n");
	    w.write("SA abs: " + absErr[3][i] + "\n");
	    w.write("\n");        
	}
	w.close();
        w = new BufferedWriter(new FileWriter(results + ".dat"));
        for (int i = 0; i < fracs; ++i)
	    {   
	    w.write(frac[i] + "");
	    for (int j = 0; j < 4; ++j)
		w.write("\t" + relErr[j][i] + "\t" + absErr[j][i]);
	    w.write("\n");
	}
        w.close();
    }

 public void varySize(String d1, String results) throws Exception // makes code simpler to avoid try/catch
    {
	int n = 10000;
	int datarepeats = 10;
	int repeats = 1;
	Random r = new Random();
	int[] ens = {1000, 10000, 100000}; //{0.05, 0.1, 0.15, 0.20, 0.25};
	int num = ens.length;
	double[][] relErr = new double[4][num];
	double[][] absErr = new double[4][num];

	for (int i = 0; i < num; ++i)
            for (int j = 0; j < 4; ++j)
                relErr[j][i] = absErr[j][i] = 0.0;

	BufferedWriter w = new BufferedWriter(new FileWriter(new File(results)));

	for (int l = 0; l < num; ++l)
        {
	    n = ens[l];
	    System.out.println("n = " + n);
	    w.write("n = " + n + "\n");
	    for (int k = 0; k < datarepeats; ++k)
	    {
		// read data from file
		BufferedReader r1 = new BufferedReader(new FileReader(d1 + n + "-" + k + ".dat"));
		n = Integer.parseInt(r1.readLine());
		double[] data = new double[n];
		for (int i = 0; i < n; i++)
		    data[i] = Double.parseDouble(r1.readLine());
		r1.close();
		
		// created sorted version of the data
		double[] dataSorted = new double[n];
		for (int i = 0; i < n; i++)
		    dataSorted[i] = data[i];
		Arrays.sort(dataSorted);
	    
		// compute the real KS statistic
		double D_real = 0.0;
		for (int i = 0; i < n; i++)
		{
		    D_real = Math.max(D_real, (i + 1.0)/n - cdfNormal(dataSorted[i], 0, 1));
		    D_real = Math.max(D_real, cdfNormal(dataSorted[i], 0, 1) - ((double)i)/n);
		}
		System.out.println("\nReal: D = " +  D_real);
		w.write("\nReal: D = " +  D_real + "\n");

		double f = 0.01;
		System.out.println("\nf = " + f);
		w.write("\nf = " + f + "\n");
		
		/** test the GK-data structure **/
		GreenwaldKhannaQuantileSummary s = new GreenwaldKhannaQuantileSummary((int)(f * n));
		for (int i = 0; i < n; ++i)
		    s.insert(data[i]);
		// compute the estimated KS statistic from GK sketch
		double D_GK = 0.0;
		for (int i = 0; i < n; i++)
		{
		    D_GK = Math.max(D_GK, (i + 1.0)/n - cdfNormal(s.getQuantile(((double)i)/n), 0, 1));
		    D_GK = Math.max(D_GK, cdfNormal(s.getQuantile(((double)i)/n), 0, 1) - ((double)i)/n);
		}
		System.out.print("GK: " + D_GK);
		System.out.print("; " + Math.abs(D_GK - D_real));
		System.out.println("; " + relError(D_GK, D_real));
		w.write("GK: " + D_GK + "; " + Math.abs(D_GK - D_real) + "; " + relError(D_GK, D_real) + "\n");
		relErr[0][l] += relError(D_GK, D_real)/datarepeats;
		absErr[0][l] += Math.abs(D_GK - D_real)/datarepeats;
		
		/** test the qDigest data structure **/                 
		QDigest s2 = new QDigest((int)(f * n)/3);  // compression factor = space/3                 
		for (int i = 0; i < n; ++i)  
		    s2.insert(data[i]);  
		// compute the estimated KS statistic from qDigest sketch      
		double D_QD = 0.0;     
		for (int i = 0; i < n; i++)     
		{ 
		    D_QD = Math.max(D_QD, (i + 1.0)/n - cdfNormal(s2.getQuantile(((double)i)/n), 0, 1));                     
		    D_QD = Math.max(D_QD, cdfNormal(s2.getQuantile(((double)i)/n), 0, 1) - ((double)i)/n);                 
		}
		System.out.print("QD: " + D_QD); 
		System.out.print("; " + Math.abs(D_QD - D_real));  
		System.out.println("; " + relError(D_QD, D_real));   
		w.write("QD: " + D_QD + "; " + Math.abs(D_QD - D_real) + "; " + relError(D_QD, D_real) + "\n");
		relErr[1][l] += relError(D_QD, D_real)/datarepeats;     
		absErr[1][l] += Math.abs(D_QD - D_real)/datarepeats; 
		
		/** test the reservoir sampling data structure **/      
		ReservoirSampler s3 = new ReservoirSampler((int)(f * n));  
		for (int i = 0; i < n; ++i)  
		    s3.insert(data[i]);                
		// compute the estimated KS statistic from reservoir sampling sketch   
		double D_RS = 0.0;    
		for (int i = 0; i < n; i++) 
		{   
		    D_RS = Math.max(D_RS, (i + 1.0)/n - cdfNormal(s3.getQuantile(((double)i)/n), 0, 1));
		    D_RS = Math.max(D_RS, cdfNormal(s3.getQuantile(((double)i)/n), 0, 1) - ((double)i)/n);
		}
		System.out.print("RS: " + D_RS);
		System.out.print("; " + Math.abs(D_RS - D_real));  
		System.out.println("; " + relError(D_RS, D_real));  
		w.write("RS: " + D_RS + "; " + Math.abs(D_RS - D_real) + "; " + relError(D_RS, D_real) + "\n");
		relErr[2][l] += relError(D_RS, D_real)/datarepeats;
		absErr[2][l] += Math.abs(D_RS - D_real)/datarepeats; 
	    
		/** Compute value from exact KS on sampled data **/
                double[] sample = s3.getSamples();
		Arrays.sort(sample);
		// compute the estimated KS statistic from sampling
                double D_SA = 0.0;
                for (int i = 0; i < sample.length; i++)
		{
		    D_SA = Math.max(D_SA, (i + 1.0)/sample.length - cdfNormal(sample[i], 0, 1));
		    D_SA = Math.max(D_SA, cdfNormal(sample[i], 0, 1) - ((double)i)/sample.length);
		}
                System.out.print("SA: " + D_SA);
                System.out.print("; " + Math.abs(D_SA - D_real));
                System.out.println("; " + relError(D_SA, D_real));
                w.write("SA: " + D_SA + "; " + Math.abs(D_SA - D_real) + "; " + relError(D_SA, D_real) + "\n");
                relErr[3][l] += relError(D_SA, D_real)/datarepeats;
                absErr[3][l] += Math.abs(D_SA - D_real)/datarepeats;
	    }
	}

	for (int i = 0; i < num; ++i)
	{
	    System.out.println("n = " + ens[i]);
	    System.out.println("GK: " + relErr[0][i] + ", " + absErr[0][i]);
	    System.out.println("QD: " + relErr[1][i] + ", " + absErr[1][i]);
	    System.out.println("RS: " + relErr[2][i] + ", " + absErr[2][i]);
	    System.out.println("SA: " + relErr[3][i] + ", " + absErr[3][i]);

	    w.write("n = " + ens[i] + "\n");
	    w.write("GK rel: " + relErr[0][i] + "\n");
	    w.write("GK abs: " + absErr[0][i] + "\n");
	    w.write("QD rel: " + relErr[1][i] + "\n");
	    w.write("QD abs: " + absErr[1][i] + "\n");
	    w.write("RS rel: " + relErr[2][i] + "\n");
	    w.write("RS abs: " + absErr[2][i] + "\n");
	    w.write("SA rel: " + relErr[3][i] + "\n");
	    w.write("SA abs: " + absErr[3][i] + "\n");
	    w.write("\n");        
	}
	w.close();
        w = new BufferedWriter(new FileWriter(results + ".dat"));
        for (int i = 0; i < num; ++i)
	{   
	    w.write(ens[i] + "");
	    for (int j = 0; j < 4; ++j)
		w.write("\t" + relErr[j][i] + "\t" + absErr[j][i]);
	    w.write("\n");
	}
        w.close();
    }

 public void varyDistance(String d1, String results) throws Exception // makes code simpler to avoid try/catch
    {
	int n = 10000;
	int datarepeats = 10;
	int repeats = 1;
	Random r = new Random();
	String[] mean = {"0.1", "0.2", "0.3", "0.4", "0.5", "0.6", "0.7", "0.8", "0.9", "1.0"}; // used strings to preserve formatting
	int means = mean.length;
	double[] real = new double[means];
	double[][] est = new double[4][means];

	for (int i = 0; i < means; ++i)
	{
	    real[i] = 0.0;
	    for (int j = 0; j < 4; ++j)
                est[j][i] = 0.0;
	}
	
	BufferedWriter w = new BufferedWriter(new FileWriter(new File(results)));

	for (int l = 0; l < means; ++l)
	{
	    System.out.println("mean = " + mean[l]);
            w.write("\nmean = " + mean[l] + "\n");


	    for (int k = 0; k < datarepeats; ++k)
	    {
		// read data from file
		BufferedReader r1 = new BufferedReader(new FileReader(d1 + mean[l] + "-1-" + k + ".dat"));
		n = Integer.parseInt(r1.readLine());
		double[] data = new double[n];
		for (int i = 0; i < n; i++)
		    data[i] = Double.parseDouble(r1.readLine());
		r1.close();
		
		// created sorted version of the data
		double[] dataSorted = new double[n];
		for (int i = 0; i < n; i++)
		    dataSorted[i] = data[i];
		Arrays.sort(dataSorted);
		
		// compute the real KS statistic
		double D_real = 0.0;
		for (int i = 0; i < n; i++)
		{
		    D_real = Math.max(D_real, (i + 1.0)/n - cdfNormal(dataSorted[i], 0, 1));
		    D_real = Math.max(D_real, cdfNormal(dataSorted[i], 0, 1) - ((double)i)/n);
		}
		System.out.println("\nReal: D = " +  D_real);
		w.write("\nReal: D = " +  D_real + "\n");
		real[l] += D_real/datarepeats;
		
		double f = 0.01;
		System.out.println("\nf = " + f);
		w.write("\nf = " + f + "\n");
		
		/** test the GK-data structure **/
		GreenwaldKhannaQuantileSummary s = new GreenwaldKhannaQuantileSummary((int)(f * n));
		for (int i = 0; i < n; ++i)
		    s.insert(data[i]);
		// compute the estimated KS statistic from GK sketch
		double D_GK = 0.0;
		for (int i = 0; i < n; i++)
	        {
		    D_GK = Math.max(D_GK, (i + 1.0)/n - cdfNormal(s.getQuantile(((double)i)/n), 0, 1));
		    D_GK = Math.max(D_GK, cdfNormal(s.getQuantile(((double)i)/n), 0, 1) - ((double)i)/n);
		}
		System.out.print("GK: " + D_GK);
		System.out.print("; " + Math.abs(D_GK - D_real));
		System.out.println("; " + relError(D_GK, D_real));
		w.write("GK: " + D_GK + "; " + Math.abs(D_GK - D_real) + "; " + relError(D_GK, D_real) + "\n");
		est[0][l] += D_GK/datarepeats;
	    	    
		/** test the qDigest data structure **/                 
		QDigest s2 = new QDigest((int)(f * n)/3);  // compression factor = space/3                 
		for (int i = 0; i < n; ++i)  
		    s2.insert(data[i]);  
		// compute the estimated KS statistic from qDigest sketch      
		double D_QD = 0.0;     
		for (int i = 0; i < n; i++)     
		{ 
		    D_QD = Math.max(D_QD, (i + 1.0)/n - cdfNormal(s2.getQuantile(((double)i)/n), 0, 1));                     
		    D_QD = Math.max(D_QD, cdfNormal(s2.getQuantile(((double)i)/n), 0, 1) - ((double)i)/n);                 
		}
		System.out.print("QD: " + D_QD); 
		System.out.print("; " + Math.abs(D_QD - D_real));  
		System.out.println("; " + relError(D_QD, D_real));   
		w.write("QD: " + D_QD + "; " + Math.abs(D_QD - D_real) + "; " + relError(D_QD, D_real) + "\n");
		est[1][l] += D_QD/datarepeats; 
				
		/** test the reservoir sampling data structure **/      
		ReservoirSampler s3 = new ReservoirSampler((int)(f * n));  
		for (int i = 0; i < n; ++i)  
		    s3.insert(data[i]);                
		// compute the estimated KS statistic from reservoir sampling sketch   
		double D_RS = 0.0;    
		for (int i = 0; i < n; i++) 
		{   
		    D_RS = Math.max(D_RS, (i + 1.0)/n - cdfNormal(s3.getQuantile(((double)i)/n), 0, 1));
		    D_RS = Math.max(D_RS, cdfNormal(s3.getQuantile(((double)i)/n), 0, 1) - ((double)i)/n);
		}
		System.out.print("RS: " + D_RS);
		System.out.print("; " + Math.abs(D_RS - D_real));  
		System.out.println("; " + relError(D_RS, D_real));  
		w.write("RS: " + D_RS + "; " + Math.abs(D_RS - D_real) + "; " + relError(D_RS, D_real) + "\n");
		est[2][l] += D_RS/datarepeats;
	    
		/** Compute value from exact KS on sampled data **/
                double[] sample = s3.getSamples();
		Arrays.sort(sample);
		// compute the estimated KS statistic from sampling
                double D_SA = 0.0;
                for (int i = 0; i < sample.length; i++)
		{
		    D_SA = Math.max(D_SA, (i + 1.0)/sample.length - cdfNormal(sample[i], 0, 1));
		    D_SA = Math.max(D_SA, cdfNormal(sample[i], 0, 1) - ((double)i)/sample.length);
		}
                System.out.print("SA: " + D_SA);
                System.out.print("; " + Math.abs(D_SA - D_real));
                System.out.println("; " + relError(D_SA, D_real));
                w.write("SA: " + D_SA + "; " + Math.abs(D_SA - D_real) + "; " + relError(D_SA, D_real) + "\n");
                est[3][l] += D_SA/datarepeats;
	    }
	}

	for (int i = 0; i < means; ++i)
	{
	    System.out.println("mean = " + mean[i]);
	    System.out.println("GK: " + est[0][i]);
	    System.out.println("QD: " + est[1][i]);
	    System.out.println("RS: " + est[2][i]);
	    System.out.println("SA: " + est[3][i]);

	    w.write("mean = " + mean[i] + "\n");
	    w.write("GK: " + est[0][i] + "\n");
	    w.write("QD: " + est[1][i] + "\n");
	    w.write("RS: " + est[2][i] + "\n");
	    w.write("SA  " + est[3][i] + "\n");
	    w.write("\n");        
	}
	w.close();
        w = new BufferedWriter(new FileWriter(results + ".dat"));
        for (int i = 0; i < means; ++i)
	{   
	    w.write(mean[i] + "\t" + real[i]);
	    for (int j = 0; j < 4; ++j)
		w.write("\t" + est[j][i]);
	    w.write("\n");
	}
        w.close();
    }

    // computes relative error of est from real
    public static double relError(double est, double real)
    {
	return Math.abs(est - real)/real;
    }
    
    // the cdf of N(mu, sigma) on value x
    public static double cdfNormal(double x, double mu, double sigma)
    {
	return 0.5 + 0.5 * erf((x - mu)/(sigma * Math.sqrt(2.0)));
    }

    // erf function taken from:
    // http://introcs.cs.princeton.edu/java/21function/ErrorFunction.java.html
    // fractional error in math formula less than 1.2 * 10 ^ -7.
    // although subject to catastrophic cancellation when z in very close to 0
    // from Chebyshev fitting formula for erf(z) from Numerical Recipes, 6.2
    public static double erf(double z) {
        double t = 1.0 / (1.0 + 0.5 * Math.abs(z));
        // use Horner's method
        double ans = 1 - t * Math.exp( -z*z   -   1.26551223 +
				       t * ( 1.00002368 +
					     t * ( 0.37409196 + 
						   t * ( 0.09678418 + 
							 t * (-0.18628806 + 
							      t * ( 0.27886807 + 
								    t * (-1.13520398 + 
									 t * ( 1.48851587 + 
									       t * (-0.82215223 + 
										    t * ( 0.17087277))))))))));
        if (z >= 0) return  ans;
        else        return -ans;
    }


}